/*
 *  cylk_regionDetection.c
 *  cylk
 *
 *  Created by Benjamin Heasly on 7/14/09.
 *  Copyright 2009 __MyCompanyName__. All rights reserved.
 *
 */

#include "cylk_regionDetection.h"

#pragma mark stack operations
int pushValueOnStack (stackItem **top, long value) {
    stackItem *newItem, *nextItem;

    // make sure that this value isn't already in the stack
    nextItem = *top;
    while (nextItem != NULL) {
        if (value == nextItem->value)
            return 0;
        else
			nextItem = nextItem->next;
    } 
	
	// got to end of stack, need to add new item
	newItem = (stackItem *) malloc(sizeof(stackItem));
	newItem->value = value;
	newItem->next = *top;
	*top = newItem;	
	return 1;		
}

long popValueOffStack (stackItem **top) {
	long value;
		
    if (*top==NULL) {
		value = (long)NULL;
    } else {
		stackItem *popped = *top;
	    *top = (popped->next);
	    value = popped->value;
        free(popped);
    }
    return(value);
}

#pragma mark region growth

int growRegionWithThresholdFromSeedXYZ(cylkImage* dataIm, cylkImage* regionIm, float regionThreshold, int x, int y, int z) {    
	float lowVoxel, highVoxel;
	short seedVoxel;
	int i, regionVoxelCount=0;
	long seedIndex;
	cylkImage thresholdIm = new_cylkImageLikeExisting(dataIm);

	fprintf(stderr, "---- grow region from voxel at x=%d, y=%d, z=%d\n", x, y, z);
	
	subscriptsToIndex(dataIm, x, y, z, &seedIndex);	
	seedVoxel = dataIm->pt[seedIndex];
	lowVoxel = seedVoxel*(1 - regionThreshold/100.0);
	highVoxel = seedVoxel*(1 + regionThreshold/100.0);
		
	// apply threshold to input data
	for(i=0; i<dataIm->nVoxels; i++)
		thresholdIm.pt[i] = ((dataIm->pt[i] >= lowVoxel) && (dataIm->pt[i] <= highVoxel));
	
	// detect region near seed voxel
	regionVoxelCount = growRegionFromSeedIndex(&thresholdIm, regionIm, seedIndex);
	
	destroy_cylkImage(&thresholdIm);
	
	fprintf(stderr, "---- ---- region has %d voxels\n", regionVoxelCount);
	return regionVoxelCount;
}

int growRegionFromSeedIndex(cylkImage* thresholdIm, cylkImage* regionIm, long seedIndex) {
    int i, numNeighbors, regionVoxelCount=0;
    long index, neighborIndexes[26];
	stackItem *stackTop = NULL;
	
	fprintf(stderr, "---- grow region from voxel at index=%d/%d\n", thresholdIm->nVoxels);
	
    if (!thresholdIm->pt[seedIndex]) {
		fprintf(stderr, "---- ---- seed voxel was not part of any image region\n");
        return -1;
    }
	
	// keep pushing neighbor voxels indexes on stack until all are redundant
	regionIm->pt[seedIndex] = 1;
	regionVoxelCount += pushValueOnStack(&stackTop, seedIndex);
	while ((long)NULL != (index = popValueOffStack(&stackTop))){
		// add neighbor voxels that are in the input image, but not yet in the grown region
		numNeighbors = indexesOfVoxelsNeighboringImageVoxel(thresholdIm, index, neighborIndexes);
		for (i=0; i<numNeighbors; i++){ 
            if (!regionIm->pt[neighborIndexes[i]] && thresholdIm->pt[neighborIndexes[i]]){
				regionIm->pt[neighborIndexes[i]] = 1;
                regionVoxelCount += pushValueOnStack(&stackTop, neighborIndexes[i]);
            }
        }
    }
    return(regionVoxelCount);
}